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Abstract 

In order to perform accurate and stable long-term numerical integration of the Einstein equations, 
several hyperbolic systems have been proposed. We here report our numerical comparisons between 
weakly hyperbolic, strongly hyperbolic, and symmetric hyperbolic systems based on Ashtekar's connec- 
tion variables. The primary advantage for using this connection formulation is that we can keep using 
the same dynamical variables for all levels of hypcrbolicity. Our numerical code demonstrates gravita- 
tional wave propagation in plane symmetric spacetimes, and we compare the accuracy of the simulation 
by monitoring the violation of the constraints. By comparing with results obtained from the weakly 
hyperbolic system, we observe the strongly and symmetric hyperbolic system show better numerical 
performance (yield less constraint violation), but not so much difference between the latter two. 
We also study asymptotically constrained systems for numerical integration of the Einstein equations, 
which are intended to be robust against perturbative errors for the free evolution of the initial data. 
First, we examine the previously proposed "A-system", which introduces artificial flows to constraint 
surfaces based on the symmetric hyperbolic formulation. We show that this system works as expected 
for the wave propagation problem in the Maxwell system and in general relativity using Ashtekar's con- 
nection formulation. Second, we propose a new mechanism to control the stability, which we call the 
"adjusted system" . This is simply obtained by adding constraint terms in the dynamical equations and 
adjusting its multipliers. We explain why a particular choice of multiplier reduces the numerical errors 
from non-positive or pure-imaginary eigenvalues of the adjusted constraint propagation equations. This 
"adjusted system" is also tested in the Maxwell system and in the Ashtekar's system. This mechanism 
affects more than the system's symmetric hypcrbolicity. 

1 Introduction 



Numerical relativity - solving the Einstein equations numerically - is now an essential field in gravity 
research. The current mainstream of numerical relativity is to demonstrate the final phase of compact 
binary objects related to gravitational wave observations and these efforts are now again shedding 
light on the mathematical structure of the Einstein equations. 

Up to a couple of years ago, the standard Arnowitt-Deser-Misner (ADM) decomposition of the 
Einstein equations was taken as the standard formulation for numerical relativists. Difficulties in accu- 
rate/stable long-term evolutions were supposed to be overcome by choosing proper gauge conditions and 
boundary conditions. Recently, however, several numerical experiments show that the standard ADM 
is not the best formulation for numerics, and finding a better formulation has become one of the main 
research topics. 

One direction in the community is to apply conformally decoupled and tracefree re-formulation of 
ADM system which were first used by Nakamura et al. [Q . Although there is an effort to show why this 
re-formulation is better than ADM, we do not yet know this method is robust for all situations. 

* Proceedings of The 10th Workshop on General Relativity and Gravitation at Osaka, Japan, September 2000. 
This report is based on our works [Q] and [^. 
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Another alternative approach to ADM is to formulate the Einstein equations to reveal hyperbolicity 
1^ . A certain kind of hyperbolicity of the dynamical equations is essential to analyze their propagation 
features mathematically, and are known to be useful in numerical approximations. Several hyperbolic 
formulations have been proposed to re-express the Einstein equations, with different levels: weakly, 
strongly and symmetric hyperbolic systems. Several numerical tests were also performed in this direction. 

The following questions, therefore, naturally present themselves (cf. 

(1) Does hyperbolicity actually contribute to the numerical accuracy /stability? 

(2) If so, which level of hyperbolic formulation is practically useful for numerical applications? (or 
does the symmetric hyperbolicity solve all the difhculties?) 

(3) Are there any other approaches to improve the accuracy/stability of the system? 

In this report, we try to answer these questions with our simple numerical experiments. Such 
comparisons are appropriate when the fundamental equations are cast in the same interface, and that is 
possible using Ashtekar's connection variables |@, More precisely, the authors' recent studies showed 
the following: 

(a) the original set of dynamical equations proposed by Ashtekar already forms a weakly hyperbolic 
system g, 

(b) by requiring additional gauge conditions or adding constraints to the dynamical equations, we can 
obtain a strongly hyperbolic system 

(c) by requiring additional gauge conditions and adding constraints to the dynamical equations, we 
can obtain a symmetric hyperbolic system |l^ , and finally 

(d) based on the above symmetric hyperbolic system, we can construct a set of dynamical systems 
which is robust against perturbative errors for constraints and reality conditions [pl]| (afca. A- 
system fl^). 

Based on the above results (a)-(c), we developed a numerical code which handles gravitational wave 
propagation in the plane symmetric spacetime. We performed the time evolutions using the above 
three levels of Ashtekar's dynamical equations together with the standard ADM equation. We compare 
these for accuracy and stability by monitoring the violation of the constraints. We also show the 
demonstrations of our A-system (above (d)), together with new proposal for controlling the stability. 

It is worth remarking that this study is the first one which shows full numerical simulations of 
Lorentzian spacetime using Ashtekar's connection variables. This research direction was suggested 
soon after Ashtekar completed his formulation, but has not yet been completed. Historically, an appli- 



cation to numerical relativity of the connection formulation was also suggested |15j using Capovilla- 
Dell-Jacobson's version of the connection variables |l6|, which produce an direct relation to Newman- 
Penrose's ^'s. However here we apply Ashtekar's original formulation, because we know how to treat its 
reality conditions in detail, and how they form hyperbolicities. We will also describe the basic numerical 
procedures in this paper. 

Due to the limited space, we have to omit details in some part in this report. More complete 
discussion can be found in our recent articles ^ . 

2 Field equations to be compared 

2.1 Hyperbolic formulations 

We begin by providing our definitions of hyperbolic systems. 

Definition 1 We assume a certain set of (complex) variables Ua {a = l,---,n) forms a first-order 
(quasi-linear) partial differential equation system, 

dtUa^f'^a{u)dlUi3 + K^{u), (1) 

where J (the characteristic matrix) and K are functions of u but do not include any derivatives of u. 
We say that the system (W is: 
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(I), weakly hyperbolic, ij all the eigenvalues of the characteristic matrix are real. 

(II). strongly (diagonalizable) hyperbolic, if the characteristic matrix is diagonalizable and has all 
real eigenvalues. 

(Ill), symmetric hyperbolic, if the characteristic matrix is a Hermitian matrix. 

We treat J'^'q as a n x n matrix when the /-index is fixed. We say A' is an eigenvalue of a when the 
characteristic equation, det {J''^ a — A'^'^q) — 0, is satisfied. The eigenvectors, p", are given by solving 

These three classes have the relation (III) £ (II) G (I). The symmetric hyperbolic system gives us 
the energy integral inequalities, which are the primary tools for studying well-posedness of the system. 

For more concrete descriptions for each systems, please refer our paper It might be worth 
remarking that the standard ADM formulation cannot be a first-order hyperbolic form, since there are 
curvature terms in the r.h.s. of dynamical equations. 



2.2 Ashtekar's dynamical equations: three levels of hyperbolicity 

The key feature of Ashtekar's formulation of general relativity is the introduction of a self-dual 
connection as one of the basic dynamical variables. The new basic variables are the densitized inverse 
triad, £^^, and the S0(3,C) self-dual connection. A";, where the indices i, j, • • ■ indicate the 3-spacetime, 
and a, 6, ■ • ■ are for SO (3) space. The total four-dimensional spacetime is described together with the 
gauge variables N,N^^Aq, which we call the densitized lapse function, shift vector and the triad lapse 
function. The system has three constraint equations, 

C^SH {^/2)e^\ElElF^^^Q, (2) 

C^f := -KM"^^^ (3) 
C^f := A^:«0, (4) 

which are called the Hamiltonian, momentum, and Gauss constraints equation, respectively. The dy- 
namical equations for a set of (i?*,^^) are 

dtK - -~iV,{e'\NEiEl) + 2V,{N^^E^)+iAleab^El, (5) 
dtAI = -le'^'.NElF^j+N^F^^+V.A^, (6) 



where F^j := 2d[iA'^^ — ie'^bcA'^Aj is the curvature 2-form. 

We have to consider the reality conditions when we use this formalism to describe the classical 
Lorentzian spacetime. Fortunately, the metric will remain on its real-valued constraint surface during 
time evolution automatically if we prepare initial data which satisfies the reality condition. More prac- 
tically, we further require that triad is real-valued. But again this reality condition appears as a gauge 
restriction on v4n |r^ , which can be imposed at every time step. In our actual simulation, we prepare 
our initial data using the standard ADM approach, so that we have no difficulties in maintaining these 
reality conditions. 

The set of dynamical equations (|^) and (^ [the original equations] does have a weakly hyperbolic 
form 1^, so that we regard the mathematical structure of the original equations as one step advanced 
from the standard ADM. Further, we can construct higher levels of hyperbolic systems by restricting 
the gauge condition and/or by adding constraint terms, C^^^jC^fj^ and Cgf^, to the original equations 
1^. We extract only the final expressions here. 

In order to obtain a symmetric hyperbolic system ^\ we add constraint terms to the right-hand-side 
of ^ and (^. The adjusted dynamical equations, 

dtK = ~iV,ie-\NEiEl) + 2V,iNi^E'J)+iA'oeab'El + KiP\i,Cf''\ (7) 
18|] p resented a symmetric hyperbolic expression in a different form. The differences between ours and theirs 



^' Iriondo et al 
are discussed in 
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where P\b = N'Sab + iN^ab^K, 
dtAt = -^e'^\NElF^J+WF^^+V,A''o + ^^2Q''^Cf + n^R^^CtfJ'. (8) 
where Q1 = e-^NEf, R,^" = ie-^Ne^^E^El 

form a symmetric hyperbohcity if we further require ki ^ K2 = K3 — 1 and the gauge conditions, 

A'i = A^N\ d^N = 0. (9) 

We remark that the adjusted coefficients, P'^ab,Qi,Ri''°', for constructing the symmetric hyperbohc 
system are uniquely determined, and there are no other additional terms (say, no C^^^,C^/^ for 
no Cq^^ for dtAf) p. The gauge conditions, (^, are consequences of the consistency with (triad) reahty 
conditions. 

We can also construct a strongly hyperbolic system by restricting to a gauge iV' ^ 0, -izNyJ^ (where 
7'' is the three- metric and we do not sum indices here) for the original equations (|^), (^. Or we can 
also construct from the adjusted equations, (0) and (||), together with the gauge condition 

A^a = A1N\ (10) 

As for the strongly hyperbolic system, we hereafter take the latter expression. 



system 


variables 


Eqs of motion 


remark 


I Ashtekar (weakly hyp.) 




( 


5), (|6 


) (original) 


"original Ashtekar" 


II Ashtekar (strongly hyp.) 




( 




) (adjusted with k = 1) 
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( 
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) (adjusted with different k) 





Table 1: List of systems that we compare in this article. 



2.3 Alternative approaches to obtain robust evolution 1: "A-system" 

Based on the symmetric hyperbolic feature of the system, Brodbeck, Frittelli, Hiibner and Reula (BFHR) 
proposed an alternative dynamical system to obtain stable evolutions, which they named "A-system" 
Ip^ . The idea of this approach is to introduce additional variables. A, which indicates the violation of 
the constraints, and to construct a symmetric hyperbolic system for both the original variables and As 
together with imposing dissipative dynamical equations for As. BFHR constructed their A-system based 
on Frittelli-Reula's symmetric hyperbolic formulation of the Einstein equations , and we ||ll| have also 
presented the similar system for the Ashtekar's connection formulation based on the above symmetric 
hyperbolic expression. Here we present our system which evolves the spacetime to the constraint surface, 
Ch ~ Cmi ~ Cca ~ as the attractor. In Jill , we also present a system which controls the perturvative 
violation of the reality condition. 

We introduce new variables (A, Ai, Aa), as they obey the dissipative evolution equations 

dtX = aiCH-f3iX, (11) 
dtXi = a2CAfi - /32 Aj, (12) 
dtXa = asCca-PsXa, (13) 

where ai ^ (allowed to be complex numbers) and f3i > (real numbers) are constants. 

If we take mL^^^ = (£'^, A, Ai, Ao) as a set of dynamical variables, then the principal part of 
(pT|)-([r^) can be written as 

dtX ^ -iaie''"'EiEl,idiA^), (14) 
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,5iEi]{diA)), 



(15) 
(16) 



The characteristic matrix of the system u'^^^ does not form a Hermitian matrix. However, if we 
modify the right-hand-side of the evolution equation of {El,Af), then the set becomes a symmet- 
ric hyperbolic system. This is done by adding 1537'' (9; Aa) to the equation of dtEl, and by adding 
a2(-e7''"£;f + eS^E''''){diX,n) to the equation of dtAf. The final principal part, 



taie^/EfEliidiX) - 
then, is written as 



9, 



' Af ^ 
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' a m 
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iaie'^/EfElj a2e{5'^E^'' - -f^'^E^) 







'E, 



El - 5\E^) 




di 



J 



A^ 
X 

Am 

Va, / 

(17) 

Clearly, the solution (£"* , A";, A, Ai, Xa) = {E^, A";, 0, 0, 0) represents the original solution of the Ashtekar 
system. If the As decay to zero after the evolution, then the solution also describes the original solution 
of the Ashtekar system in that stage. Since the dynamical system of Ua^^\ (p^, constitutes a symmetric 
hyperbolic form, the solutions to the A-system are unique. BFHR showed analytically that such a decay 
of As can be seen for As sufhciently close to zero with a choice of appropriate combination of as and [3s, 



and that statement can be also applied to our system. Therefore, the dynamical system, (17), is useful 



for stabilizing numerical simulations from the point that it recovers the constraint surface automatically. 



2.4 Alternative approaches to obtain robust evolution 2: "adjusted system" 

We also try to compare a set of evolution system, which we propose as "adjusted-system" . The essential 
procedures are to add constraint terms in the right-hand-side of the dynamical equations with multipliers, 
and to choose the multipliers so as to decrease the violation of constraint equation. This second step will 
be explained by obtaining non-positive (or non-zero) eigenvalues of the adjusted constraint propagation 
equations. We remark that this eigenvalue criterion is also the core part of the theoretical support of 
the above A-system. 

The fundamental equations that wc will demonstrate in this report, are the same with (0) and (|8|), but 
here the real-valued constant multipliers ks are not necessary equals to unity. We set k = ki — K2 — H'i 
for simplicity. Apparently the set of (j^) and (^) becomes the original weakly hyperbolic system if k = 0, 
becomes the symmetric hyperbolic system if k = 1 and N = const.. The set remains strongly hyperbolic 
systems for other choices of n except k = 1/2 which only forms weakly hyperbolic system. 



3 Comparing numerical performance 
3.1 Model and Numerical method 

The model we present here is gravitational wave propagation in a planar spacetime under periodic 
boundary condition. We perform a full numerical simulation using Ashtekar's variables. We prepare two 
-|— mode strong pulse waves initially by solving the ADM Hamiltonian constraint equation, using York- 
O'Murchadha's conformal approach. Then we transform the initial Cauchy data (3-metric and extrinsic 
curvature) into the connection variables, (-E^,^°), and evolve them using the dynamical equations. For 
the presentation in this article, we apply the geodesic slicing condition (ADM lapse TV = 1 or densitized 
lapse = 1, with zero shift and zero triad lapse). We have used both the Brailovskaya integration 

scheme, which is a second order predictor-corrector method, and the so-called iterative Crank-Nicholson 
integration scheme for numerical time evolutions. The details of the numerical method are described in 
PI, where we also described how our code shows second order convergence behaviour. 
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More specifically, we set our initial guess 3-metric as 



7y 






l + /f(e-(^--^)' +e-(^'+'^)') |, (18) 



in the periodically bounded region x = [—5, +5]. Here K and L are constants and we set K = 0.3 and 
L = 2.5 for the plots. 

In order to show the expected "stabilization behaviour" clearly, we artificially add an error in the 
middle of the time evolution. We added an artificial inconsistent rescaling once at time t = 6 for the Ay 

■y -^yK 



component as ^ „4^(1 + error). 



3.2 Differences between three levels of hyperbolicity 

We have performed comparisons of stability and/or accuracy between weakly and strongly hyperbolic 
systems, and between weakly and symmetric hyperbolic systems[0. (We can not compare strongly and 
symmetric hyperbolic systems directly, because these two requires different gauge conditions.) 

We omit figures in this report, but one can see a part of results in Fig.|l| and Fig.||. We may 
conclude that higher level hyperbolic system gives us slightly accurate evolutions. However, if we evaluate 
the magnitude of L2 norms, then we also conclude that there is no measurable differences between 
strongly and symmetric hyperbolicities. This last fact will be supported more affirmatively in the next 
experiments. 



3.3 Demonstrating "A-system" 

Next, we show a result of the "A-system" Fig.|l| (a) shows how the violation of the Hamiltonian 
constraint equation, Ch, become worse depending on the term error. The oscillation of the L2 norm Ch 
in the figure due to the pulse waves collide periodically in the numerical region. We, then, fix the error 
term as a 20% spike, and try to evolve the same data in different equations of motion, i.e., the original 
Ashtekar's equation [solid line in Fig.|l| (b)], strongly hyperbolic version of Ashtekar's equation (dotted 
line) and the above A-system equation (other lines) with different /3s but the same a. As we expected, 
all the A-system cases result in reducing the Hamiltonian constraint errors. 



3.4 Demonstrating "adjusted system" 

We here examine how the adjusted multipliers contribute to the system's stability |Q. In Fig.^, we show 
the results of this experiment. We plot the violation of the constraint equations both Ch and Cmx- An 
artificial error term was added in the same way as above. The solid line is the case of k = 0, that is 
the case of "no adjusted" original Ashtekar equation (weakly hyperbolic system). The dotted line is 
for K — I, equivalent to the symmetric hyperbolic system. We see other line (k = 2.0) shows better 
performance than the symmetric hyperbolic case. 



4 Why adjusted equations have better performance? 

The remaining question is: why we can get the better performance by adding constraint terms in the 
dynamical equations? The added terms are basically error terms during the evolution for its original 
dynamical equations. Nevertheless, these terms improve the accuracy of the evolution. This question 
applies to both higher level of hyperbolicities and also for our proposing "adjusted systems" . Here we 
introduce our plausible explanation schematically. The detail explanations and numerical experiments 
are in 1^. 

Suppose we have constraint equations, Ci « 0, C2 « 0, • • •, in a system. We, normally, monitor the 
error of the evolution by evaluating these constraint equations on the each constant-time hypersurface. 
Such monitoring, on the other hand, can be performed also by checking the evolution equations of 
the constraint, which we denote constraint propagation equations (cf. |2l| ). If this set of constraint 
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Figure 1: Demonstration of tlie A-system in tire Aslitekar equation. We plot tlie violation of the constraint (L2 norm 
of the Hamiltonian constraint equation, Ch) for the cases of plane wave propagation under the periodic boundary. 
To see the effect more clearly, we added artificial error at t = 6. Fig. (a) shows how the system goes bad depending 
on the amplitude of artificial error. The error was of the form Ay -^yi^ + error). All the lines are of the evolution 
of Ashtekar's original equation (no A-system). Fig. (b) shows the effect of A-system. All the lines are 20% error 
amplitude, but shows the difference of evolution equations. The solid line is for Ashtekar's original equation (the 
same as in Fig. (a)), the dotted line is for the strongly hyperbolic Ashtekar's equation. Other lines are of A-systems, 
which produces better performance than that of the strongly hyperbolic system. 



propagation equations could be written in a first order form, then we may predict the evolution behavior 
by its characteristic matrix, M, which is expressed by 




where ~ means we are extracting only the characteristic part. 

The idea here is to estimate the eigenvalues of the characteristic matrix, M, after we took the Fourier 
transformation on the both sides of (p^. 

1. Clearly, if all the real part of the eigenvalues are negative, then all constraints decays to zero along 
to the system's evolution. 

2. Alternatively, if all the eigenvalues are pure-imaginary, then all constraints evolve similar to wave 
equations, and diffusion due to differenciation may preserve from its unstable evolution. 

These guidelines are supported by von Neumann's stability analysis. For the adjusted Ashtekar equa- 
tions, the eigenvalues of the amplicication matrix in von Neumann's method for the set of {Ch, Cmi, Cgo) 
are {1,1 — v^nsinO ± sin 0| (multiplicity 2 each),l — i/^(l — 2K)^sin^6' ± iv\l — 2K||sin6'|}, where 
V = At/ Ax, 9 is the parameter for the frequency in grid spacing domain, and n is the adjusted param- 
eter (§2.4). We observe that von Neumann's necessary condition for the stability can be obtained in 
larger non-zero k. 

In 1^, we show that such a case can be obtained by adding 'adjusted terms' both for Ashtekar's 
and Maxwell's systems. There we also show examples of unstable evolution by choosing adjusted terms 
which produce positive eigenvalues of M . 

In this point, we can say that adjusted terms are responsible for obtaining the stable and/or accurate 
evolution system, and this is a way to control the stability of simulation, which effects more than the 
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Figure 2: Demonstration of the adjusted system in the Ashtekar equation. An artificial error term was added 
at i = 6, in the form of Ay -^^(1 + error), where error is +20% as before. Fig. (a) and (b) are L2 norm of 
the Hamiltonian constraint equation, C//, and momentum constraint equation, Cmx, respectively. The solid line is 
the case of k = 0, that is the case of "no adjusted" original Ashtekar equation (weakly hyperbolic system). The 
dotted line is for k — 1, equivalent to the symmetric hyperbolic system. We see other line (k = 2.0) shows better 
performance than the symmetric hyperbolic case. 



system's hyperbolicity. 



5 Concluding remarks 

We presented numerical comparisons on the accuracy/stability of the systems, mainly between three 
mathematical levels of hyperbolicity: weakly hyperbolic, strongly hyperbolic, and symmetric hyperbolic 
systems. We apply Ashtekar's connection formulation, because this is the only known system in which we 
can compare three hyperbolic levels with the same interface. We also tested two another approaches with 
a purpose of constructing an "asymptotically constrained" system, "A-system" and "adjusted system" , 
which can be robust against perturbative errors for the free evolution of the initial data. 

We may conclude that higher levels of hyperbolic formulations help the numerics more, though its 
differences are small Two other approaches work as desired. These show quite good performance 
than those of the original and its symmetric hyperbolic form. This indicates, in turn, that the symmetric 
hyperbolic system is not always the best for controlling accuracy or stability. 

The reason why higher hyperbolicity shows better stability may be explained by the eigenvalue 
analysis of the propagation equation of the constraints. Our two guidelines, negative-real eigenvalues or 
pure-imaginary eigenvalues, give us clear indications whether the constraints decay (i.e. stable system) 
or not for perturbative errors, at least to all our numerical experiments. 

We think these results open a new direction to numerical relativists for future treatment of the 
Einstein equations. 

To conclude, we are glad to announce that Ashtekar's connection variables have finally been applied 
in numerical simulations. This new approach, we hope, will contribute to understanding further of 
gravitational physics, and will open a new window for peeling off interesting non-linear natures together 
with a step to numerical treatment of quantum gravity. 
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